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Abstract 

m 

^ ' We compare the eigenvalue spectra of the Dirac operator from a simulation with two mass de- 

generate dynamical chirally improved fermions with Random Matrix Theory. Comparisons with 
distribution of k— th eigenvalues (k — 1,2) in fixed topological sectors (y — 0,1) are carried out 
using the Kolmogorov-Smirnov test. The eigenvalue distributions are well described by the RMT 
predictions. The match allows us to read off the quark condensate in the chiral limit. Correcting for 
finite size and renormalization we obtain a mean value of -(276(11)(16) MeV) 3 in the MS scheme. 

O 

1 Introduction 

Ph" _ - . 

The light quark condensate £ = (uu) « (dd) is a measure of chiral symmetry breaking. Most prominently 
it is the proportionality factor between the renormalized quark mass and the experimentally measurable 
(/ttM^) 2 (the Gell-Mann-Oakes-Renner relation [T]) in the chiral limit. In the effective description of low 
energy chiral symmetry breaking (chiral perturbation theory [2]) it is a fundamental externally supplied 
parameter. Other features, like the microscopic eigenvalue distribution of the Dirac operator, follow 
universal laws related just to the general symmetry structure - there the condensate is the only parameter 
which is to be provided. Only full QCD calculations include the necessary dynamics to determine this 
parameter ab initio. 

QCD breaks chiral symmetry spontaneously, the small quark masses provide an additional explicit break- 
ing. The condensate is affected by both, but as far as is known, only weakly by the explicit breaking, i.e., 
it has a distinctive non-zero value in the chiral limit. Both, the physical quark mass and the condensate 
are renormalization prescription dependent and have to be given in some scheme. In the continuum MS 
scheme (at scale 2 GeV) the condensate has been determined from QCD sum rules to values between 
E 1 / 3 = —250 ... — 270 MeV [3]. A recent large-N expansion gives similar values, not very sensitive to the 
number of flavors 0]. 

Lattice methods to compute £ include measurements of ratios of correlation functions (for recent deter- 
minations in the framework of the Dirac operator used here, the Chirally Improved (CI) operator, see 
[S]). These, however, need lattices of sufficient size to identify the asymptotic behavior of the correlators 
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and from these the renormalizcd quark mass. Random matrix theory, on the other hand, allows the 
determination on small lattices. 

Older global averages, using two dynamical flavors of staggered, Wilson and clover improvement type 
of fermions generated values like E 1 / 3 = —280(13) MeV for the chiral condensate [B]. More recent 
determinations using parameters of chiral Lagrangians obtained from lattice QCD [7] yielded values like 
£i/3 _ _259(27) MeV in the MS at (2 GeV) for 2+1 flavors of dynamical staggered quarks (cf. that 
reference also for a comparison with other authors' values). Domain wall fermions have also been used to 
measure the condensate for two dynamical flavors [5], giving values of the unrenormalized E 1 / 3 between 
-223 MeV and -261 MeV. 

Studies with dynamical overlap fermions on small (10 4 ) lattices, using RMT distribution like we do here, 
report values of E 1 / 3 = -269(9) MeV for N f = 1 and E 1 / 3 = -282(10) MeV for N f = 2 QD]. The 
latter values are without a finite volume (lattice shape) correction factor, which would lower the value to 
-269 MeV0 

A recent study with staggered sea quarks and overlap valence quarks obtain an unrenormalized value of 
—291(5) MeV for E 1 / 3 [TT], with estimated renormalization and finite volume corrections to be 13 - 18%. 
Our own studies with 2 dynamical chirally improved fermions, based on the GMOR relation, indicated a 
value of E 1 / 3 = -288(8) MeV [H]. 

Here we study spectra of the CI Dirac operators, obtained in our dynamical simulation, comparing 
the resulting density distributions with the RMT parameterizations in order to obtain values for the 
condensate including finite size corrections. The approach is similar to that in [TU], the difference being 
mainly a different Dirac operator and larger lattices, here at four different combinations of quark mass 
and lattice spacing. 

2 Technical background 

2.1 Simulation with CI fermions 

The CI Dirac operator is an approximate Ginsparg- Wilson operator [12! ; it is a truncated series solution 
[T3"| to the GW relation. CI fermions have been already extensively tested in quenched calculations (see, 
e.g., Ref. [14 ). There it was found that one can reach pion masses below 300 MeV without running 
into the problem of exceptional configurations (spurious zero modes) . On quenched configurations pion 
masses down to 280 MeV could be obtained on lattices of size 16 3 x 32 (lattice spacing 0.148 fm). 

In recent works [TBI \T7\ [T8] we have studied the CI fermions in a dynamical simulation of QCD with two 
light flavors. All technicalities are discussed in Ref. [T7]. We use the Luscher-Weisz gauge action [T5] and 
stout smearing 20J of the gauge fields as part of the Dirac operator definition. The Hybrid Monte Carlo 
method was implemented to deal with the dynamics of the fermions. The lattices were (up to now) of 
moderate size: 8 3 x 16 and 12 3 x 24, with lattice spacings between 0.11 and 0.14 fm. Table [T] summarizes 
the simulation parameters of the runs discussed here (see, however, Ref. 117] for a more complete list and 
details of the simulation). 

In the table we also give the number of configurations of the run sequences where the low lying spectrum 
of the Dirac operator has been extracted. Typically this was done every 5th configuration in order to 
reduce autocorrelation effects. Simulations with dynamical fermions are very costly, even more so for the 
GW-type fermionic actions. For the moment we therefore have to rely on comparatively small samples 
and lattices. 

2.2 Random Matrix Theory 

Random Matrix Theory (RMT) was introduced into physics by Wigner [3T] who used it to describe the 
apparently universal distribution of nuclear level spacings. The connection between RMT and the low 

1 After completion of this paper we learned of recent work [32| presenting a two-flavor simulation with overlap quarks in 
the e-regime. The result quoted there is -251(7)(11) MeV. 
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run 


L 3 


x T 


Pi 


am 


HMC 

time 


as[fm] 


a niAwi 


conf 


P 


Zs,q 


a 


12 3 


x 24 


5.2 


0.02 


463 


0.115(6) 


0.025 


73 


1.21 


1.07 


b 


12 3 


x 24 


5.2 


0.03 


363 


0.125(6) 


0.037 


52 


1.18 


1.10 


c 


12 3 


x 24 


5.3 


0.04 


438 


0.120(4) 


0.037 


55 


1.19 


1.08 


d 


12 3 


x 24 


5.3 


0.05 


302 


0.129(1) 


0.050 


40 


1.17 


1.10 



Table 1: Parameters for the simulations; the first column denotes the run, for later reference. The gauge 
coupling is (3\, the bare quark mass parameter am, HMC-time denotes the length of the run (number 
of trajectories), and the lattice spacing has been determined via the Sommer parameter. In the last 
two columns we give the finite volume correction parameter p and the renormalization factor for the 
condensate obtained in a quenched setting [H], as discussed in Sect. [3] 

lying spectrum of the Dirac operator was made in |22[ I23j . Since then several lattice studies (mostly for 
the quenched situation) [24] have confirmed the predictions of RMT. (See [25, 2(3 [27] for recent reviews 
in that context.) 

The present understanding is that the generating function for universal features of the low lying Dirac 
spectrum is determined by chiral symmetry. The microscopic spectrum thus follows the spectral properties 
of a random matrix with symmetry like the Dirac operator, i.e., the so-called chiral Gaussian Unitary 
Ensemble (chGUE). This agreement is expected to hold in a regime where the pion wave length is larger 
than the volume considered, such that only pions dominate the dynamics and all other heavy degrees 
of freedom become irrelevant. In chiral perturbation theory this is the so-called the e-regime 28 . Our 
pion mass is clearly too large for that with m„ » 3/L; however, the smallest values of Im(A) leading to 
RMT distributions have values such that ■\/2|Im(A)I]|// 7r m l/L and thus we may hope that RMT is 
applicable. 

The QCD parameters, quark mass and condensate, enter the RMT distributions in form of dimensionless 
scaling variables /i = Y.Vm and £ = T,VX, where A denotes eigenvalues of the matrix (or the imaginary 
part of the Dirac operator eigenvalues, respectively), and V is the volume. The procedure for calculating 
the exact RMT distributions Pfc[^](C) of the fcth largest eigenvalue for a given number of flavors and 
in a fixed topological sector v have been worked out and checked in [29]. In this way RMT allows one 
to determine E from comparison of numerically computed eigenvalue distributions with the universal 
distributions given in terms of the scaling variables. Since our operator is not exactly chiral, we define 
our masses through the Axial Ward Identity (the so-called AWI-mass) to take into account possible 
additive renormalization. 

Even though the predictions of RMT are true only in the infinite volume limit with m tending to zero 
keeping fj, fixed, the distributions are approached already for comparatively small systems (for recent 
applications cf. PUHH]). The obtained value of E, although corresponding to the parameter of infinitely 
sized random matrices, is not yet the condensate in a lattice simulation. There the geometry (hypercubic, 
different space- and time extent) provides extra correction factors which may be computed from sub- 
leading chiral perturbation theory. One also has to consider scale dependent renormalization factors. We 
discuss both in Sect. [3] 

As emphasized, e.g., in |10| , the bias introduced by histogramming the data may be avoided using the 
Kolmogorov-Smirnov (KS) test (cf. [30] for details). To apply the KS test one first creates a cumulative 
distribution function of both the data as well as the probability function from which it is supposed to 
be drawn and compares the two. We will use it to test whether the RMT distribution, for fixed v and 
fermion mass, for some value of E describes the distribution of the individual eigenvalues of the CI Dirac 
operator. Given a data set, we find the best fit to the RMT distribution with the single parameter "S 1 / 3 " 
(which minimizes the distance D between the data cumulative distribution and the RMT expectation) 
and quote the probability we obtain for that value of E 1 / 3 . 
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Figure 1: Normalized histograms and cumulative distributions for run (a). The corresponding prob- 
abilities (Q K s(k, v)) are : Q KS (1, 0)=0.50, QksQ-, 1)=0.99, Qks(2, 0)=0.90, Qks(2, 1)=0.96. The 
abscissa C — EVA, where A denotes the imaginary part of the corresponding eigenvalues. 
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Figure 2: Like Fig. [TJ but now for run (d): Qk S (1, 0)=0.98, Qks(1, 1)=0.99, Qks(2, 0)=0.74, 
Q ks (2, 1)=0.99. 

3 Results: Condensate 

For exact GW-operators the eigenvalues are exactly on the GW-circle, which may be projected to the 
imaginary axis for further analysis, as done in [10 . In our case the eigenvalues typically are close to, but 
not exactly on the circle. For a sample of the localization see Ref. [TH]. We therefore use the distribution 
of the values Im(A) for comparing with the RMT cumulative distribution. The approximation introduces 
a systematic error. This adds to the distortion effects due to lattice shape and finite volume and thus 
we cannot expect our data to approach perfect RMT distributions. We expect, however, that the central 
scaling properties of the peak position survive in good approximation. 

For the mass parameter of the RMT distributions we use the AWI-mass as determined in [T7| . Only real 
eigenmodes carry chirality and these are used to determine the topological sector v = n_ — n+, where 
n± denotes the number of real modes with right or left chirality. In most (97%) configurations all the 
real modes have the same chirality. 

In Figs.Q]andE]we show some of the histograms and cumulative distributions together with the theoretical 
cumulative distributions for the optimal value of E. 

The results for all cumulants that have been analyzed are given in Table [U By E 1 / 3 we denote the 
values extracted from the distributions according to the RMT formulas. The given values maximize the 
KS-probability for the distribution. The discrepancy between the sum of the number of values entering 
the histograms and the total number of configurations analyzed is due to the configuration in higher 
topological sectors. 
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run 


k v 


# values 


-S x /3[MeV] 


Qks 


-Dobs 


(a) 


1 


19 


325 


( 25) 


0.50 


0.18 




1 1 


27 


280 


(5) 


0.99 


0.08 




2 


19 


288 


(3) 


0.93 


0.12 




2 1 


27 


284 


(3) 


0.97 


0.09 




averaj 


;e: 


285 


(3)(21) 






(b) 


1 


8 


279 


( io) 


1.00 


0.08 




1 1 


20 


338 


(11) 


1.00 


0.08 




2 


8 


290 


(6) 


0.77 


0.22 




2 1 


20 


313 


(9) 


0.78 


0.14 




averaj 


;e: 


301 


(20) (23) 






(c) 


1 


16 


291 


( 20 ) 


0.62 


0.18 




1 1 


21 


305 


(9) 


0.87 


0.13 




2 


1G 


289 


(6) 


0.99 


0.11 




2 1 


21 


301 


(7) 


0.29 


0.21 




averaj 


;e: 


294 


(7)(7) 






(d) 


1 


12 


259 


(4) 


0.98 


0.13 




1 1 


12 


300 


(6) 


1.00 


0.09 




2 


12 


286 


(4) 


0.74 


0.19 




2 1 


12 


294 


(2) 


1.00 


0.09 




averaj 


;e: 


288 


(13)(16) 







Table 2: Results for the optimal value of E obtained from fitting the individual cumulants to the RMT 
distributions, adjusting the parameter E 1 / 3 . The runs are ordered as in Table [TJ In the table k denotes 
the k th zero in the topological sector v (0 or 1). The number of eigenvalues entering the individual 
distributions is given (# values), as well as the KS-probability and the best value of the distance D bs- 

The statistical errors for E 1 / 3 given in the table have been calculated by using statistical bootstrap on 
the data. From each set of data for a given run, k, and v we produce several sets which are reanalyzed to 
give the spread in the resulting condensate. For each sector we have just few eigenvalues and therefore 
the statistical fluctuation is significant and may be even underestimated. 

We choose to fit the distributions for each k, v set individually. A joint fit would somewhat obscure the 
qualitative differences which we did observe. Naively, RMT is expected to give a unique E for all these 
samples for large enough statistics. However, as mentioned, we cannot expect exact RMT distribution 
shapes even for large enough statistics. In view of these limitations the observed differences of values 
from the individual fits given in Table [5] are not surprising. From the variation of the values we may 
derive an estimate for a systematic error (which does not cover all possible sources like scaling with the 
lattices spacing and properties of the Dirac operator). 

The sensitivity of the KS tests for a cumulative distribution function P(x) is not independent of x. In 
fact the KS test tends to be most sensitive around the median value P(x) = 0.5 and less sensitive at the 
extreme ends of the distribution where P(x) is near or 1. The result is that while KS tests are usually 
good for finding shifts in probability distribution, they are not so good at finding spreads which generally 
affect the tails more than the median. Physically, in our case, this means that a high confidence level in 
the KS tests indicates that the position of the maximum in the eigenvalue distribution is well predicted 
by RMT, but does not tell us much about the width of the distributions. 

In order to give an overall estimate we average our central values for E 1 / 3 for the distributions for each 
run, weighted by the KS-probability and the inverse statistical error squared, leading to the numbers 
given in the summary lines of Table [2l A precise definition of the statistical error is hardly possible 
due to the interplay of different probability factors. We estimate the statistical error from the variance 
around the central value, again respecting the individual statistical errors and the KS-probabilitics. The 
first error following the mean value denotes the statistical one and the second provides an estimate of the 
systematic error as discussed. 

For a check of the quality of our error estimates we also compute the best fit values of E 1 / 3 for k=3, v=Q 
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and see if they lie within the error bars of the mean values for each lattice as given in Table [21 For the 
runs a, 6, c, d we obtain for E 1 / 3 , the values 290(4), 290(2), 288(13), 281(4) with Q KS of 0.98, 0.76, 0.13, 
0.48 respectively. All these values are consistent within statistical errors with the mean values obtained 
from the fits to k=l ,2 and v=0, 1. 

Leading order chiral perturbation theory provides corrections to the condensate E depending of the 
lattice shape and the volume. The geometry (hypercubic, different space- and time extent) provides extra 
correction factors which have been discussed in sub-leading chiral perturbation theory. The condensate 
E(V) as obtained from the finite volume simulation has to be divided by a correction factor p given by 

[sii m 

N 2 f -i a 

E(oo) = E(7)/p with p = l+ I— „ 2 ■ (1) 

N f f4 L 

Here L = V 1 ' d , f n is the pion decay constant, and 3 depends upon the lattice geometry. The procedure 
for calculating this coefficient in general is outlined in [31] . Here we only quote the result relevant for our 
geometry, i.e., for a lattice twice as long in time direction than in the other three: 3 = 0.0836011. 

The values of f„ for runs (a)-(d) have been measured [17] . The extrapolation to the chiral limit is 
compatible with the physical value 93 MeV, within a 10% error margin. Since the correction factor 
applies to the chiral limit, we can estimate it using the physical value for the decay constant and the 
lattice spacings as determined in [T7J . We give the factors in Table Q] 

One also has to consider scale dependent renormalization factors Z$ — 1/Z m in order to relate to, e.g., 
the MS scheme values, = Zg E. For the quenched case these have been determined in [TS] for various 
lattice spacings, in the chiral limit, leading to values ranging from 1.13 (at a = 0.148 fm) to 0.96 (at 
a = 0.078 fm). We quote the factors Zg,q obtained from an interpolation of these quenched values in 
Table [TJ For the lattice spacings studied here they are all close to 1.1. We have not determined them 
for the dynamical simulation (where an extrapolation to the chiral limit would not be very stable) but 
expect similar values. 

Table [5] gives the uncorrected mean values for the condensate. In order to compare with the continuum 
MS values we therefore multiply the E values of the table with (Zs/p). For the runs (a)-(d) this results 
in final values 

EW(a) = (-274(3) (20) MeV) 3 , 
S (r) (t) = (-293(20) (23) MeV) 3 , 
sW(c) = (-285(7) (7) MeV) 3 , 

E w (d) = (-282(13)(16) MeV) 3 . (2) 

The weighted average for the mean, with a simple average for the errors gives an overall conservative 
estimate of E — (— 276(11)(16) MeV) 3 . This is compatible with determinations by other groups and also 
with our own determination based on the GMOR relation. 



4 Discussion and Conclusion 

In this article we have extracted the quark condensate by comparing the low lying eigenvalue spectra 
of CI-Dirac operator with Random Matrix Theory. The advantage of this method is that it gives the 
condensate directly in the chiral limit already on moderate sized lattices. By maximizing the probability 
in the Kolmogorov-Smirnov tests we obtained the chiral condensate for each of the k(— 1, 2) eigenvalue 
distributions in fixed topological sectors v(= 0, 1). Our final estimate after correcting for renormalization 
and finite volume effects is (ipip) = (— 276(11)(16) MeV) 3 which is consistent with determinations by other 
authors. Random Matrix Theory seems to predict the peaks of the individual distributions quite well. 
The widths of the distributions on the other hand are not so well matched. 
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